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, Abstract. We simulate convection near the solar surface, where the 

^Sj ' continuum optical depth is of order unity. Hence, to determine the ra- 

diative heating and cooling in the energy conservation equation, we must 
solve the radiative transfer equation (instead of using the diffusion or 
optically thin cooling approximations). A method efficient enough to cal- 
culate the radiation for thousands of time steps is needed. We assume 
If^ . LTE and a non-gray opacity grouped into 4 bins according to strength. 

On I We perform a formal solution of the Feautrier equation along a vertical 

. and four straight, slanted, rays (at four azimuthal angles which are ro- 

Q I tated 15 deg. every time step). We present details of our method. We 

also give some results: comparing simulated and observed line profiles 
for the Sun, showing the importance of 3D transfer for the structure of 
the mean atmosphere and the eigenfrequencies of p- modes, illustrating 
O I Stokes profiles for micropores, and analyzing the effect of radiation on 

-4—' ' p-mode asymmetries. 



1. Introduction 

Radiative energy exchange is critical in determining the structure of the upper 
convection zone. Near the surface of the Sun, the energy flow changes from al- 
most exclusively convective below the surface to radiative above the surface. The 
interaction between convection and radiation near the surface determines what 
we observe, drives the convection, and generates the magnetic and non-magnetic 
activity which heats the chromosphere and corona. Hence, the interaction be- 
tween convection and radiation has significant impact on both the dynamics 
of convection and our diagnostics. Escaping radiation cools the plasma that 
reaches the surface, which produces the low entropy plasma whose buoyancy 
work drives the convection. Radiation transport determines (with convection 
and waves) the mean atmospheric structure. Radiation provides the diagnostics 
that we use to determine the velocity, density, temperature and magnetic field 
on the Sun. Radiation transport modifies the observable p-mode temperature 
fluctuations so as to reverse the asymmetry of the intensity vs. velocity spectral 
peaks. 
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Since the top of the convection zone occurs near the level where the contin- 
uum optical depth is one, neither the optically thin nor the diffusion approxima- 
tions give reasonable results. We need to solve the 3D, LTE, non-gray radiation 
transfer in our models (Nordlund 1982). 

2. Radiative Heating/Cooling 

Radiation enters the calculation through the radiative heating/cooling term in 
the energy equation. 



+ V • Ve = V • P + Qrad + Qdissipation , (1) 



de ^ P 
+ V • Ve = 

ot p 

where the radiative heating/cooling is given by 



Qrsd = ATr I^Kx{Jx-Sx)dX. (2) 

We calculate Jx — Sx hy a formal solution of the Feautrier equations assuming 
LTE so the source function is the Planck function. 



d^P 



drx 



^ = Px-Bx 



where P\ is the Feautrier mean intensity 



Px = \[h{^)+ix{-m 



3. Opacities 

Continuous opacity is calculated explicitly using the Uppsala package (Gustafs- 

son et al. 1975), for the entire temperature-density range in the simulation and 
stored in a table. To account for UV opacities not in the original package, the 
UV continuous opacity is enhanced below 500 nm by a function of A, according 
to the prescription of Magain (1983). Line opacities for the mean ID atmosphere 
state are calculated using the Uppsala ODF tables (Gustafsson et al. 1975). For 
a grid of temperatures of the mean atmosphere state, the Rosseland mean with 
and without lines is calculated. The ratio of full to continuous Rosseland opac- 
ities is then used to enhance the continuous Rosseland opacity at each density 
in the table for the given temperature. In addition, for small optical depths, 
the continuum Rosseland mean opacity is corrected by the ratio of the intensity 
weighted opacity for the ID plane parallel average atmosphere to the continuum 
Rosseland opacity. Thus the opacity stored in the table is 

K = Xft^KcR , (3) 

where the correction factor Xn{T) is 

= e-^°-- #^ + f 1 - e-'^^^A . (4) 

(«)cR ^ ^ ('«)cR 
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The continuum Rosseland opacity is 

1 ^ j{i=l) 



where j{i = 1) is the set of all those wavelengths belonging to the continuum. 
The intensity mean is 

j 

To summarize: the opacity table is calculated from the Rosseland mean 
opacity, based on the continuum opacity alone, {k)cr- This is then corrected 
by the ratio (Eqn 4), to account for lines and the transition to optically 
thin radiation, but this ratio is evaluated only for the p/T-track of the mean 
stratification of the simulation (not for the whole table). The correction factor, 
Xk,, is then regarded as a function of temperature alone and used to correct 
the continuum Rosseland opacities for all densities in the table at each table 
temperature. 
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Figure 1. Opacity as a function of wavelength (left) and grouped into 

four bins (right). Each bin has the average value of the opacity in that 
bin and the weight equal to the sum of the weights of the wavelengths 
belonging to that bin. 

The radiation calculation is greatly sped up by drastically reducing the 
number of wavelengths at which we solve the Feautrier equations. We make 
the simplifying approximation of grouping the opacity at each wavelength into 
4 bins according to its magnitude (Nordlund 1982, Skartlien 2000). Suppose the 
opacity is as shown in the toy example (Fig. 1). Each wavelength is assigned 
a bin according to the magnitude of its opacity. In this toy example, those 
wavelengths with opacity between 2.0 and 2.5 are assigned to the first opacity 
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bin with average opacity 2.3 and considered to represent the continuum (Fig 1). 
Those wavelengths with opacity between 2.5 and 3.0 are assigned to bin two 
with average opacity 2.7 corresponding to weak Unes. Those wavelengths with 
opacity between 3.0 and 3.5 are assigned to bin three with average opacity 3.2 
corresponding to intermediate strength lines. Those wavelengths with opacity 
greater than 3.5 are assigned to bin 4 with average opacity 3.8 corresponding to 
strong lines. In practice, wc take the opacity bins to have a ratio of a factor of 
10 in magnitude (Nordlund 1982), 



The weight of each opacity bin is the sum of the weights of each wavelength in 
that bin, 

Wi = ^j(i)Wxj , 

where is the non-contiguous set of those wavelengths Xj in bin i. The con- 
sequence of this grouping of opacities is that all wavelengths in a given bin have 
optical depth unity at approximately the same geometrical depth so that inte- 
grals over optical depth commute with the sum over wavelengths. Specifically, 
the radiation heating/cooling is 



Qrad = j^Kx {Jx - Bx) dX 



= E E '^A, ( Ja, - 5a J u;a, 

= EE'^A.-^^r,. {Bx,)wx, 

= J2^ri{Bi)Wi 

i 

= E ^-^i ~ ^ 



Wi 



where the operator L-,-^ is 
and Bi is the weighted sum of the Planck functions at the wavelengths in bin i, 



Jx-Bx = Lr,{Bx) = /' ^e--/^ r dte-'I^Bx{t) - Bx , 
Jo u Jo 



= E-^A,^^^A,■ M 



This results in a tremendous saving in computational time, since the Feautrier 

equations need only be solved for a few opacity groups. Wc choose to use only 
4 groups corresponding to the continuum, weak, medium and strong lines. 
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Figure 2. The Feautrier equation is solved along one vertical and 4 
inclined rays through each point on the surface. Values of the opacity 
and source function are interpolated to the intersection of the ray with 
each horizontal plane. The azimuthal angle (p of the inclined rays is 
rotated 15 degrees each time step in order to sample the surrounding 
medium. 

A major approximation we make is to assume that bin membership of each 
wavelength is independent of temperature. This is reasonable to the extent 
that the ratio of line to continuum opacity is relatively independent of the tem- 
perature, which holds for lines of atomic species that are also major electron 
contributors to the dominant continuous opacity. However, it is not a good 
approximation for the numerous iron lines that produce significant line blocking 
in the photosphere. 

4. Solution of Feautrier Equations 

We solve the Feautrier equations along 5 rays - one vertical ray and four slanted 
rays (at an angle cos 6 = /j, = 1/3), through each grid point on the surface, using 
long characteristics (Fig. 2). The vertical ray has weight 1/4 and the slanted 
rays together weight 3/4. The four slanted rays are rotated by 15 degrees each 
time step in order effectively probe all the surrounding medium. Each 0, cf) 
combination is treated in turn. For each, there is a collection of parallel rays 
labeled by the surface grid point through which they pass. For the slanted 
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rays, the opacity and source function are interpolated to the location of the 
intersections of the rays with each horizontal plane. The optical depth along 
each ray at each horizontal plane is then calculated by integrating along the ray 
from the surface inward. Since the opacity groups are fixed multiples (10) of 
each other, the optical depth needs to be calculated for the continuum only and 
is then scaled for the other opacity groups. 

To obtain accurate integrated intensities the continuum opacity average is 
renormalized slightly, with a factor (close to unity) obtained by requiring that 
the emergent intensity, summed over the four bins, should match the wavelength 
integrated monochromatic intensity, computed in a two-dimensional strip from 
the model. A posteriori checks, where similar comparisons are made for 3-D 
snapshots, have been used to verify the procedure. 

It is crucial that the thin surface thermal boundary layer be resolved in 
order for the correct radiative losses to be obtained. Therefore, in low resolu- 
tion simulations we interpolate from the hydrodynamic computational grid to a 
finer grid in the vertical direction for solving the radiation transport. In high 
resolution simulations this is not necessary because the thermal boundary layer 
is already resolved on the dynamic grid. 

At large optical depths, Px — > -B^, which leads to roundoff errors in the 
heating qx = P\ — B\. We therefore solve the Feautrier equation directly for 
q\ = P\ — Bx, by rewriting it as 



dr'x 



qx 



fBx 
drl 



(5) 



This ensures that the radiative heating/cooling goes to the correct asymptotic 
solution at large optical depth. The boundary conditions are 



drx 



0, 



qx + Bx 
(r- 



dBx 

drx 
oo) 



(r-0) 



The finite difference form of the Feautrier equation is. 



f I 2 




qj 



V^j+i - 




- 0?-l 






2 




1 "0+1 




1 




2 



Jj+l-Tj Tj+l -Tj-l 

This has the tri-diagonal form 



1 + 



1 1 

+ 



Tj - Tj-l Tj + i - Tjl Tj + l - Tj-l 



1 1 
+ 



Jj - "0-1 "O+i - "O / '0+1 - '0-1 



AjQj-i + BjQj + CjQj+l = Dj 
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There is a roundoff error problem at small optical depth, because the 1 in 
Bj becomes so small compared to that it is lost in the noise and the 

solution becomes indeterminate. As a result, this formulation can only be used 
down to optical depths the order of the square root of the machine word length. 
The solution is to store the sum of the tridiagonal matrix elements, rather than 
Bj, in the diagonal slot (Nordlund 1982, Rybicki Sz Hummer 1991), initializing 
it with 

Dj==A,+Bj + C, = -l 

The tridiagonal system is solved by the standard elimination procedure, except 
that the diagonal element Bj is recovered just before the back-substitution step, 
where it may be obtained from 

since at that point Aj = from the forward elimination. 

With these substitutions, the solution is accurate for arbitrarily small opti- 
cal depths. 

Total radiative heating/cooling at each location where a ray crosses a hori- 
zontal coordinate surface is accumulated as the sum over opacity bins {i = 1,4), 
and angles of the Feautrier qi{9,(j)), 

Qrad = 47r ^ ^ Kiqi{n)WiWn . (6) 

a i 

The radiative heating/cooling is finally interpolated from the slanted grid back 
to the Cartesian computational grid and added into the expression for the time 
derivative of the internal energy and the variables are advanced in time. 

5. Interaction of Radiation and Convection 

The interaction between radiative and convective energy transport has a pro- 
found effect on the structure, dynamics and diagnostics of the Sun. Escaping 
radiation cools the plasma that reaches the solar surface and produces the low 
entropy plasma whose buoyancy work drives the convection. Radiative, con- 
vective and wave energy transport determines the mean structure of the solar 
atmosphere. Radiation emerging from the atmosphere provides the diagnostics 
that we use to determine the velocity, density, temperature and magnetic field 
on the Sun. Radiation transport modifies the observable p-mode temperature 
fluctuations so as to reverse the asymmetry of the intensity vs. velocity spectral 
peaks. 

5.1. Convection Driver 

Solar convection is driven primarily from a thin surface thermal boundary layer 
by radiative cooling. Fluid that approaches the surface loses its thermal and 
ionization energy to escaping radiation (Fig 3). This reduces its entropy, so 
gravity pulls it down to form the cores of the downflows. Large entropy fluctu- 
ations occur only in these downflows, so that most of the buoyancy work that 
drives the convection occurs in these downflows. Thus the primary driver of 
solar convection is radiative cooling from the surface (Stein &; Nordlund 1998). 
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Figure 3. Fluid that approaches the surface is cooled by radiating 
away its thermal and ionization energy. This decrease its entropy and 
it is pulled down by gravity to form the cores of downfiow plumes. The 
largest entropy fluctuations and most of the buoyancy work that drives 
the convection occurs in the downflows. 



5.2. Atmospheric Structure 

The average structure of the solar atmosphere is significantly different from that 
of ID models that reproduce the observed solar luminosity and radius. First, due 
to 3D radiative transfer effects the average temperature of a 3D model is higher 
than a ID model with the same effective temperature. The great temperature 
sensitivity of the H~ opacity makes the optical depth of hot gas very large, so it 
is not observed (Fig 4). Hence, the actual average gas temperature is higher than 
the observed radiation temperature. A hotter atmosphere has a larger pressure 
scale height and is more extended. In the solar case this gives an extra 75 km 
extension to the atmosphere. Second, turbulent pressure gives extra support to 
the atmosphere, which extends it an additional 75 km. (Fig 5) (Rosenthal et al. 
1999). 

5.3. Atmospheric Diagnostics 

Emerging continuum and line radiation provides us with the diagnostic tools 
to determine the properties of the solar atmosphere. Emergent radiation from 
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Figure 4. Radiation temperature in the IR at 1.6 fim vs. gas temper- 
ature at < TiR >= 1. Hot gas is not observed. The large temperature 
sensitivity of the H~ opacity leads to large optical depths where the 
gas is hot and hides it from view. Thus the average gas temperature 
is higher than the observed radiation temperature, so the 3D atmo- 
sphere is more extended than a ID atmosphere with the same effective 
temperature. 



1 .000 0.999 0.998 




-0.5 0.0 0.5 1.0 1.5 2.0 

z[Mm] 



Figure 5. Atmospheric stratification is extended 150 km, half by 
3D radiative transfer effects (Fig 4) and half by turbulent pressure 
support. 

the solar surface can be resolved spatially and temporally. In the continuum 

radiation we see the granulation. The spatial spectrum (Fig 6) gives the power 
in the different scales of motion. The intensity distribution tells us about the 
temperature contrast on the unit optical depth surface. 

Line radiation gives us much more information. Spectral lines of heavy 
elements, whose thermal Doppler widths are small compared to typical photo- 
spheric velocities, provide direct diagnostics of velocity and temperature fluc- 
tuations in the photosphere. Non-spatially resolved properties such as average 
line profiles are useful because they bypass the difficulties associated with at- 
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Figure 6. Granulation size spectrum and emergent intensity distri- 
bution compared to observed values. 




Doppler velocity [km/s] Doppler velocity [km/s] 

Figure 7. Average line profile for Fe A 541.4 line, with temperature 
structure from the simulation, with no velocities (left) and with the 
full velocity field (right). 



mospheric seeing and instrumental resolution. Their shape critically depends 
on the solar velocity field - the temperature structure alone gives lines that 
are too narrow and deep (Fig 7 (left)). Including the convective overshoot ve- 
locities in the photosphere gives excellent agreement with the observed profiles 
(Fig 7 (right)). This observed average profile is the result of spatial and tem- 
poral averaging of lines from different locations with different temperatures and 
line of sight velocities and very different shifts, widths and shapes (Fig 8). The 
average line profile results entirely from convection induced temperature and 
velocity fluctuations, without any need for micro- or macro- turbulence or extra 
damping (Stein & Nordlund 2000). 

The surface magnetic field can also be determined by analyzing stokes com- 
ponent line profiles. An example of the stokes profiles along a slit across a 
micropore in the simulation is shown in Fig. 9. 

5.4. P-Modes 

The solar p-modes are excited by the work of turbulent pressure and non- 
adiabatic gas pressure (entropy) fiuctuations near the top of the convection 
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Figure 8. The average line profile (heavy line and Fig 7) is the com- 
bination of profiles from different locations with different line shifts and 
widths. The effects that are modeled by micro- and macro- turbulence 
in ID are all due to the overshooting convection velocities. 




Figure 9. Stokes profiles of Fcl 630.1 nm (right) along a slit across a 
micropore, whose intensity image is shown on the left. 

zone (Fig 10). The turbulent pressure is due to the convective motions. The 
entropy fluctuations result from the instantaneous local imbalance between the 
convective heating and radiative cooling. P-mode excitation decreases at low 
frequencies because of the increasing mode mass and decreasing mode compres- 
sion as the frequency decreases. Excitation decreases at high frequencies because 
the pressure fluctuations produced by the convection fall off with increasing fre- 
quency. The modes are excited fairly close to the solar surface - the higher the 
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frequency the closer to the surface is their excitation (Fig 10) (Stein & Nordlund 
2001). 




Figure 10. Solar p-modes are excited by the work of turbulent pres- 
sure and non-adiabatic gas pressure fluctuations. The excitation is 
confined closer to the surface as the frequency increases. 

The velocity and intensity spectra of the p-modes have asymmetric peaks. 
The velocity has more power on the low frequency side of the peak and the 
intensity has more power on the high frequency side of the peak. The velocity 
and intensity measured at a fixed geometrical depth have the same asymmetry, 
which depends on the locations of the excitation source and the observations 
(Kumar & Basu 1999, Georgobiani et al. 2000). The intensity asymmetry 
gets reversed from that of the temperature by radiation effects. The dominant 
H~ opacity is very temperature sensitive. The oscillation induced temperature 
fluctuations are larger on the low frequency side of the resonance and produce 
opacity changes that vary the location of radiation emission (r = 1) in a way 
that reduces the magnitude of the temperature fluctuations more on the low 
frequency side of the resonance, and thus reverses the asymmetry of the observed 
temperature fluctuations (Figs 11, 12). 

6. Discussion 

In order to make realistic numerical simulations it is necessary to include the 
relevant realistic physics. In the case of the solar atmosphere this means includ- 
ing an appropriately realistic treatment of the radiative losses. Many techniques 
are available to solve for radiation transport. The problem is to devise a method 
that is efficient enough to be applied at each time step in a dynamic calculation. 

We have described a method developed by Nordlund that relies on reducing 
the number of frequencies that need to be treated to a bare minimum of 4, by 
using a multi-group binning of the wavelengths according to their opacity. In 
addition a minimal number of rays are employed, but they are rotated to cover 
the computational domain. Because we are only interested in the photosphere, 
we also can assume LTE. Some of the approximations we use, such as assuming 
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Figure 11. Velocity (solid) and temperature (dotted) spectra at fixed 
geometrical depth < r >= 1 (left) and optical depth r = 1 (right), for 
the first non-radial fundamental mode (corresponding to solar mode 
with i = 740). Radiation transfer reduces the observed temperature 
fluctuations by over an order of magnitude and reverses the asymmetry 
of the temperature profile. 
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Figure 12. Opacity (left) and height of optical depth unity (right). 
Larger temperature fluctuations (dotted) at frequencies below the 
mode peak generate larger opacity fluctuations, which in turn produce 
larger fluctuations in the height of the r = 1 surface. This decreases 
the temperature fluctuations more below the peak than above it and 
reverses the asymmetry of the intensity with respect to the velocity. 
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that all the opacities have the same depth dependence, can be easily improved. 
It was instituted in a time when computer memory was small and the table 
size needed to be minimized. It would be nice to replace the approximation 
of multi-group opacities with an opacity sampling method, so as to include the 
effect of Doppler shifts on the radiation absorption, but the cost will be many 
more frequencies that need to be solved. For the Sun we are lucky that the 
radiation heating/cooling time scale is of the same order as the dynamical time 
scale, so that an explicit solution of the radiation transport is possible. For 
giant stars, for instance, where the dynamical time scale is much longer than 
the Sun's, explicit solution of the radiation transport will reduce the possible 
time step size. 

For other investigations, other physics will be needed. To include the chro- 
mosphere it is necessary to include at least the dominant non-LTE effects of 
scattering and hydrogen ionization. This is the next project we are working on 
in collaboration with Mats Carlsson, Viggo Hansteen and Andrew McMurray of 
Oslo University. 

Acknowledgments. The work of RFS was supported in part by NASA 

grant NAG 5 9563, NSF grants AST 9819799 and ATM 9988111. AN was 
supported in part by the Danish Research Foundation, through its establishment 
of the Theoretical Astrophysics Center. The calculations were performed at 
the National Center for Supercomputer Applications, which is supported by 
the National Science Foundation, at Michigan State University and at UNI»C, 
Denmark. This valuable support is greatly appreciated. 

References 

Georgobiani, D. C, Kosovichev, A. G., Nigam, R., Nordlund, A., & Stein, R. F. 

2000, ApJ, 530, L139 
Gustafsson, B., Bell, R. Eriksson, K. & Nordlund, A. 1975, A&A, 42, 407 
Kumar, P. k Basu, S. 1999, ApJ, 519, 389 
Magain, R 1983, A&A, 122, 225 
Nordlund, A, 1982, A&A, 107, 1 

Rosenthal, C. S., Cliristensen-Dalsgaard, J., Nordlund, A., Stein, R. F., & 

Trampedach, R. 1999, A&A, 351, 689 
Rybicki, G. B. & Hummer, D. G. 1991, AkA, 245, 171 
Skartlien, R. 2000, ApJ, 536, 465 
Stein, R. F. k Nordlund, A. 1998, ApJ, 499, 914 
Stein, R. F. k Nordlund, A. 2000, Sol. Phys., 192, 91 
Stein, R. F. k Nordlund, A. 2001, ApJ, 546, 585 



